function [W,normParam]=ecogWienerMIMORegressionTrain(ecog,Y,intervalLength,intervalOffset,regressionParam)
%[W]=ecogWienerMIMORegressionTrain(ecog,Y,intervalLength,intervalOffset,regressionParam) Calclulate a MIMO Wiener regression 
%
% Purpose: Calculate the Multiple Input Multiple Output (MIMO) Wiener regression 
% matrix W that does Y=X*W, i.e. predicts Y (the movement or stimulus feature)
% from X (the ecog data). Ridge regression is used to estimate W.
%
% Y is regressed with multiple shifted version of X. Therefore, W extends in
% time and establishes a certain temporal relation between X and Y that can 
% be interpreted in terms of causality (see below). Heres X causal on Y means
% that effects in X have an influence on what happens later in time in Y. 
% 
% INPUT:
% ecog:     An ecog structure with fields:
%           data
%           sampDur
% Y:        The matrix of time series of external variables (measured movements,
%           or stimulus parameters). Variables change along the first
%           dimension and time increases along the second.
% intervalLength:
%           The length of the intervall of the regression filter 
%           (in units of ecog.sampDur). The number of estimated regression parameters 
%           per external variable is intervalLengthInSamples*numberOfChannelsInX 
% intervalOffset:
%           Time offset between the brain data and external variables:
%           Four cases can be distinguished:
% 
%           intervalOffset<=-intervalLength: 
% returns a regression marix W with purely  
% causal effects of X on Y. If both absolute values are equal returns the  
% W with the shortest lag with purely causal effects of X on Y. Decreasing  
% the offset further (its a negative value!) produces X on Y causal Ws with  
% longer lags
% 
%           0>=intervalOffset>-intervalLength: 
% returns a matrix W producing mixed X 
% and Y on X causal effects. 
%
%           intervalOffset>=1: 
%produces the filter with the shortest lag with purely 
% causal effects of Y on X. Further increasing the offset (its a positive 
% value!) produces Y on X causal Ws with longer lags.
%
%           intervalOffset=0 and intervalLength=1
%           standard regression
%
% regressionParam (optional): struct with fields
%           method - 0 for ridge regression (default) 
%                    1 for 'LARS-EN' 
%                    2 for LARS
%                    3 for LASSO (see
%                    http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=3897)
%                    normalization required!!
%           lambda (default = 0);
%           normalize - 0 (default) no normalization
%                       1 zscore X 
%                       2 zscore X and Y
%           (Attention: if normalizing, keep normParam for prediction!)
%           
%
% OUTPUT:
% W:        The matrix W to predict Y from X in the equation Yhat=X*W.
%           The first enttry in W is an offset and must be discarde for
%           prediction.
% normParam: The zscore param for Y and X 
%           If regressionParam.normalize=true, use the normParam as input
%           for ecogWienerMIMORegressionPredict
%
% Requirements: 
% Currently the ststistics toolbox for ridge regression.
% LARS-EN,LARS,LASSO implememtation:
% http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=3897
% for elastic net computation

% 091129 JR wrote it based on Brian's code
% 100607 CR: improved chopping (use most possible timesteps and remove bug
% for case 1), zscoring added, elastic net calculation prepared
% 111111 CR: create sub functions ecogPrepareRegression,
% ecogTrainRegression, ecogPredictRegression
% TODO: Potential new code ridge regression code is appended at the end of the function.  


%% Input check
if nargin <5,
    regressionParam=struct;
end
intervalLengthSamp=round(intervalLength/ecog.sampDur); % ms to samples
intervalOffsetSamp=round(intervalOffset/ecog.sampDur);

[X,Y]=ecogPrepareRegression(ecog.data,Y,intervalLengthSamp,intervalOffsetSamp);

[W,normParam] = ecogTrainRegression(X,Y,regressionParam);
